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1 . INTRODUCTION 


The study of the ionosphere using the Faraday rotation effect has been 
undertaken recently by means of rocket, satellite, and moon echo experiments. 
Different approximations have been used by different authors, resulting in 
methods with different degrees of complexity, and it is possible to say that 
the more accurate a method is, the more difficult its application becomes. 

However, the use of modern high-speed digital computers offers the 
possibility of using more complex methods in the solution of this problem. 

The program described in this report was written for the ILLIAC, the 
digital computer of the University of Illinois. Only the general features 
common to most digital computers will be mentioned. 

This program was prepared having in mind the analysis of the Faraday 
rotation effect, as recorded from artificial satellites. It is intended 
to be as general as possible in the conditions imposed on the assumed 
propagating medium: specifically there are no restrictions on the models of 
the electron density distribution and the earth’s magnetic field. as. long as 
the ray theory is valid. 

The program will trace separately the ordinary and the extraordinary 
mode, and it will find the virtual phase path length of a ray of each mode 
between the transmitter (satellite) and a receiver (station) . 

, The difference between respective phase path-lengths is related to the 
Faraday rotation through a constant which depends on the frequency. 
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2, GENERAL THEORY OF FARADAY ROTATION 


As early as 1845., Faraday had discovered that the plane of polarization 
of light is rotated through a certain angle when it traverses certain sub- 
stances in the presence of a longitudinal .magnetic field.. The same effect . 
has been observed on radio waves traversing the ionosphere or any ionized 
medium with a magnetic field present, a so-called magneto-ionic medium, 
even in the .case when- the magnetic field is not exactly longitudinal. 

The reason , for this effect is the existence of two characteristic modes 
of propagation. Each of these modes has its own polarization, phase velocity 
and group velocity. The polarization and the velocities are completely defined . 
by the characteristics of; the medium, so that an originally linearly-polarized 
wave will be resolved into the two characteristic modes. The total polar- 
ization will, therefore, be different at different points in the medium. 

If a uniform magneto-ionic medium with a longitudinal field is considered, 
the characteristic modes are circularly polarized with opposite senses of . 
rotation and slightly different phase velocities. Therefore, the resultant 
polarization is linear with its plane of polarization determined by the phase 
relation .between the modes. The relation between these phases changes 
continuously along the ray, resulting, in a rotation of the plane of 
polarization.. This rotation is given by (Pedersen, 1927 1 ) . 





( 2 . 1 ) 


where ft = rotation of the plane of polarization in radians along a distance ^ 

i - distance in the phase propagation direction 

Vp = phase velocity of the ordinary mode 

v (x) _ phase velocity of the extraordinary mode 
P 

w = angular frequency of the wave 

In the ionospheric study using the Faraday effect, none of the conditions 
mentioned above are present. The radio wave is neither propagating under a 
strict longitudinal condition nor is the medium homogeneous. Two generalizations 
are necessary to the original concept of Faraday rotation expressed by (2.1) , , 
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When the longitudinal condition is not met, the resultant polarization 
is generally elliptical. Its shape, as well as its orientation, changes 
continuously along the propagation direction. IJt is not exactly a rotation 
and a careful interpretation shbuld.be made. However, ( 2 . 1 ) gives the correct 
relation between the phases of both modes. From this relation the resultant 
polarization can be deduced. In the present work it is desirable to 
generalize the concept of rotation rather than to use the quasi-long- 
itudinal and quasi-transverse: approximations. Following this point of view, 

( 2 . 1 ) is valid for the general case, and even for the transverse, case, if 

it is decided in the former that the transition from circular to linear 

TT 

polarization is equivalent to a rotation of — in £2 (Daniels and Bauer, 

2 2 
1959 ) . 


When the medium is not uniform, (2.1) will be written in differential 

dl dl 

form. For reasons of convenience 7 — r- and 7 — r will be substituted by 

/ \ / N (o) (x) y 

JrT ,(o) J Jm (x) ^ . V V 

dT and dT respectively. p p 


d ft = | (dT (o) - dT (x) ) 


( 2 . 2 ) 


By integrating this expression between two points T and R: 


co 

^TR = 2 




(2.3) 


Where is the generalized Faraday rotation between the points T and R and 

the integrals, are the times required for phase propagation of the respective 
modes between the same two points. This new generalization is made by defining 
the relation between the phases of both modes. The need of this definition 
arises- from the fact that the characteristic modes follow a different path 
in a varying medium. However, this expression for Faraday rotation is not 
arbitrary since it reflects the result that would be obtained if an observation 
point moved from T to R, investigating the resultant polarization at each point. 

In the expression (2.3), it has been assumed that both modes are defined 
at any point in the integral path as if the medium were uniform in the 
surroundings of such a point. The assumption is possible only if the change 
of the physical conditions of the medium is negligible in a distance of a 
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few wavelengths. Such a medium is called a "slowly varying medium". It is 
defined in relation to a given frequency. Fortunately, the ionosphere is 
a '"slowly varying" medium with respect to the frequencies used in the 
Faraday rotation experiments. Furthermore, the integrals can be evaluated 

independently of one another when both modes propagate independently. 

3 

^Booker, (1936 ) has shown that in a "slowly varying medium, " such as the 
ionosphere, the ordinary and extraordinary mode propagate independently of 
one another and the ray theory applies to each mode. The polarization 
ellipse of each mode changes slowly along its path, in accordance with 
local conditions (electron density, earth* s magnetic field and direction of 
propagation) . 

Two conditions are listed under which the independence of the modes 
breaks down and under which coupling between both modes should be considered: 

a) In which the refractive indices (1^°^ and |j/ X ^ approach zero, that is, 

ox » 

when the critical plasma frequencies w , w are nearly equal, to the frequency 
considered. In this case the physical' conditions vary slowly but the wave- 
length becomes very large so that the "slowly varying" condition is violated. 

A wave solution should be sought. Since satellites use frequencies well above 
the maximum plasma frequency in the ionosphere the refractive index is well 
above zero. This condition does not occur. 

b) Near the "edges" of the ionosphere, where the electron density is 
low. In this situation the independence of the modes is not preserved. The 
coupled waves of the other mode, produced over a long part of the ray, have a 
cumulative effect on the considered ray. The only result is that the 
characteristic polarization ellipse stops changing at a certain "limiting 
polarization height". This height is above the lower edge of the ionosphere. 

This height is dependent on the coupling coefficient, which in turn depends on 

4 

electron density, magnetic* field and collision frequency of electrons (Budden, 1952 ) . 

Neither the path of the ray nor the group and phase velocities of either 
mode are affected because of the coupling in this region. Therefore, the 
total Faraday rotation defined in (2 C 3) remains unchanged by this effect. 

However, if the polarization of ~the total wave at the station is examined, 
this effect should be taken into, account since the total wave is formed by the 
recombination of the ordinary and extraordinary modes as they reach the station, 
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and the polarization of each mode remains unchanged below the mentioned 
" limiting polarisation height 1 . 1 . 

The integrals which appear in Equation (2.3) between two points arise 
from integration along the respective group trajectories or ray paths. 

Figure 1 illustrates the relation among the differential of length ds along 
the ray path, the differential time dT, and the phase velocity corresponding 
to any of the modes. 

The phase velocities are given by 


v (o) 
p 




(2.4) 


where c = velocity of light in free space 

= refractive indices given by the Appleton-Hartree formula 
(Equation 3.1) and calculated in the wave-normal direction 
By observing Figure 1 and using relations (2.4) it is possible to write: 


dT 

dT 


(o) 

1 

(o) 

(o) 

~ c 


cos a 

(x) 

1 

(x) 

(x) 

” c 


cos a 


ds 


(o) 


ds 


(x) 


(2.5) 


Equations (2.5) show the convenience of using a vector |i of magnitude |i 
and with its direction in the direction of the wave normal. The final form 
of Equation (2.3) will be: 





R 

r 

^ T 




J 




( 2 . 6 ) 


Some simplified forms of Equation (.2.2) are in common use in the literature. 
These are given below: 

(a) Both modes are assumed to follow the same path as determined by the 

isotropic refractive index. In this case Equation (2.6) reduces to (Little, 

5 

and Lawrence,. .1960 .) . * . 





) ds 


(2.7) 
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Relation between d 5 dT and v 

s' p 


Figure l c 




(b) Both modes are assumed to follow a straight line connecting the 

transmitter and receiver. In this case we can further reduce Equation (2.7) 

0 

and write it in the following form (Browne, Evans, and Hargreaves, 1959 ) 



NMdh 


N = electron density, 

7 

M = H sec i (Yeh and Gonzales, 1960 ) 
L 

dh = differential height 


( 2 . 


where 



/ 
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3. THE REFRACTIVE INDEX 


g 

The Appleton-Hartree (1932) formula for the refractive index of a 
magneto-ionic medium is well known in its physical meaning; however, its 
use is sometimes obscure. Before proceeding to a presentation of the mathe- 
matical method, this topic will be briefly examined. If .the collision 
frequency is disregarded^ 




(o) 


X 

If 


1 


X 


1 


Y' 2 sin 2 ® / Y^sin 4 © 
2<1 ' !0 "V 4(1-X) 2 


+ Y 2 cos 2 0 


(3.1) 


(1 


(o) 




(x) 


0 


refractive ^.ndex for the ordinary mode, corresponding to the 
upper sign. 

refractive index for the extraordinary mode, corresponding to 
the lower sign. 

angle between wave normal and magnetic field. 


X = 


2 

e 


2 

4t r me 


o 



Y = 


e|i o 
27 T m 


H 

f 


(3.2) 


e = charge of an electron (in coulombs) 

m = mass of an electron (in kilograms) 

6 = dielectric constant of free space (MKS) 

o 

ll = magnetic permeability of free space (MKS) 
o 

f = frequency of wave under consideration (cycles per second) 

N = electron density (electrons per cubic meter) 

H = intensity of magnetic field (amperes per meter) 

The upper index or for the ordinary or extraordinary mode will 

be omitted on (i and on the other quantities, when the equations could refer 
to either of the modes. 
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From the foregoing, itiiis possible to see that for a given frequency, 
f, the only quantities not constant are N, H, and ©. In the steady state, 

N will be a scalar function of point (x,y,z, or r, 0, <P) . H will be a 
vector function of the coordinates, that is, the magnitude and the 
direction of H will be defined in terms of x, y,z, or r, 0, according to 
the coordinate system chosen. Therefore, it is possible to write (3.1) in 
functional form in the following way: 


n = |i(x,y,z, 0) 


for rectangular coordinates 


or 

11 = |i(r, 0,^,0) 


for spherical coordinates 


(3.3) 


The angle © is the angle between two directions, that of the magnetic 
field and that of the wave normal. The direction of the magnetic field is 
determined by the variables x, y, z (or r, 0 , . To determine the direction 

of the wave normal two angles are needed, such as azimuth and altitude, 
referred to a local orthogonal system centered at the point x, y, z(or r, 

e, <p). 

It is preferable, instead of two angles, to use the three that the wave 

normal will form with each of the local orthogonal axes: i , i , i (or i , 

x * y * z v’ 

ig, i with an additional condition.: 


2 2 2 
cos i + cos i + cos i = 1 
x y z 


2 2 2 
cos. i^ + cos ig + cos i^ = 1 


(3.4) 


In the method that follows, a vector of components v X , v^, v z (or v r , 

Q (p 

v , v ) referred to a local orthogonal system centered at the point x, y, 
z, (or r, 0, ty) has been chosen to determine the direction of the wave 
normal, subject to the following condition: 


2 

x 


v 


+ V" 


2 

z 

+ V 



r 2 0 2 2 


(3.5) 
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In this way (3.3) can be written in functional form as: 

IX = H-(x,y,z,v X ,v y ,v Z ) 

r 0 (Q 
fJ. = H'(r,0,<P v ^ v , v 1 ) 


(3.6) 


From condition (3.5) , the following relations are verified: 


= H* cos i 


= M* cos i 


= H* cos i 


= M* cos i 


or 


- H* cos i 


0 


.<p 


= |x cos i. 


(3.7) 


L <P 

For the purposes of the present report and from the last relations, the 
v T s can be considered as the component of a vector H* of magnitude H-. How- 
ever, it should be remembered that their principal function is to define a 
direction. 

A surface on which the refractive index in a given direction is constant 

\ 

is defined as a "stratification surface". For example: 


H-(x, y, z, v X , v y , v Z ) = const 

v = const 
x 

v = const 

y 

v - const 
z 


(3.8) 


In Section 5 it will be shown that Snell's law is verified with respect 
to these stratification surfaces in the differential ray tracing Equations (4.1). 

The surfaces (3.8) are not the same as: 

r 0 (0 

(->-(r,0,^,v ,v ,v ) = const 


v = const 

0 

v = const 

<p 

v = const 


(3.9) 


r 0 

Since the same set v , v , ^v define a different direction in each point of 

space, the partial derivatives-^) x v z > if ] x v z )/ ) xv 

x / y> z ; v , v ,v y / x,z,v ,v,v z / x, y, v , v. 
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will represent the components of a vector perpendicular to the stratification 
surfaces with similar characteristics to those of a gradient. These partial 
derivatives would represent the variation of ^ in space when H- is considered 
in the same direction at every point. The corresponding relations in spherical 
coordinates do not have the same characteristics because a change of the 
coordinates 0 and <P with the v* s constant implies a change of direction in space 
for H*. For example, if a uniform medium is considered, i.e., a medium with 
constant electron density and constant magnetic field (in magnitude and 
direction) through all space, then: 


v 

r 

X 


x y z ~&— 
y,z,v ,v ,v y 


V 

w 

z 


) 


X y Z 

x iy, v , v , v 


= 0 


(3.10) 


But if a spherical coordinate system is chosen: 


S 

H- 

cT 

r 




e cp 


= 0 ; 




r,(p } v X 


e <p 






) a r e (p 
Jr,Q,v ,v ,v 


t 0 (3.11) 


A final point to consider is the fact that the stratification surfaces 
do not correspond to surfaces of constant electron density, unless the 
magnetic field is constant over the surface. 
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4„ MATHEMATICAL METHOD 


4^1' Ray Differential Equations 

The ray tracing equations used are those derived from Hamilton’s method 
and expressed in general, rectangular and spherical coordinates by Haselfrove 
(1955) 9 . 

The following set of equations for spherical coordinates was chosen; 


dr 

eft 



av 


d0 

dt 



- H- 


3^ . 
av 


d (p 
dt 


r sin 0 H" 


(0 

(v r - 


V- 


0H- . 
av r 


dy r _ _1 3H- 
dt~ ~ V- 9r 


0 d0 

v — + 
dt 


sin 0 v 


<P &P 

dt 


(4.1) 


dv 

dt 


.i { i 

r L H* 


an 

W 


~ v 


9 dr 
dt 


+ r cos 


(p d <p 


dt 


dv 

dt 


<P 


r sin 


j i 

1 v- 


an- 

wp 


- sin 


<P dr 
dt 


v 


- r cos 


<P d0 

7 dt 


where; H* = refractive index in the wave normal direction 

r, 0 <P = spherical coordinates referred to a system with its origin 

at the center of the Earth. The z axis passes through the 

north pole. The Greenwich meridian is in the xz plane. 

T* 0 (fi ■ — * 

v ' 9 v 9 v ss components of along the local orthogonal system associated 
with the point (r^ 6^ <P) 
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t = virtual phase path length. This variable is related to the time 

of phase propagation along the ray path by the velocity of light c. 

The integration of these equations gives the trajectory of the "group" 
path or "ray" corresponding to the "o" or "x” mode according to the chosen 
value of the refractive index. The final value of t is the phase length 
along the ray. 

The partial derivatives — ^ } have to be transformed into 

3v- 


9v 


3v 


a more workable form since their meaning is somewhat obscure. 

2 * 0 (D 

The functional relation of with respect to v ; v and v can be 
expressed ; * V. 


n = n [r, e, % O] 


© = © ( V y, v^, v^) 


= cos 


r 0 tp 

_ V COS *] + V COS + V COS 'Hn 

r u u 


/ 


r 2 e 2 

v + v + v 


(4.2) 


where cos *1^, cos 'Hq, cos ^ = direction cosines of the earth* s magnetic field 
referred to the local orthogonal system, associated with the point (r, 0, <f) . 

i 

< ^-Kr^n^pia'tXbns' (4.2): 

^ ^© 


9v r 


9v r 


9© 


9p _ 

d v 90 3 V ® 


(4.3) 


9p _ 9jj 90 
9 ^ ~ 9 0 9 ^ 

Performing operations (4.3) on expressions (4.2) and substituting into 
the first three of equations (4.1), 


dr 1 r ** _ . r /Ci 

~ ~ [ v - D (v cos ^ - fi cos *1 ) ] 

ii r 


(4.4) 
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d9 1,0 n , © 0 „ 

[v - D (v cos - p. cos Tig) ] 


dt 2 
r(i 


d <P 

dt 


(p (p r\ 

[v - D (v cos v? - ji cos r l < ^) ] 


r sin 9 p. 


(4.4) 


where D is defined as; 

1 ik 


D = 


(X sin © §0 


2 ft 

Y cos ^ 


4 4 

I y sin © 2 Q 

r _r - ft v 


2_U i_ 

2(1-X) 


+ for (o) 
- for (x) 


Y cos 


4(1 -x) " 


(4.5) 


The set of differential equations used to obtain* the ray are set (4.4) 

and the last three equations of (4*1) „ 

Condition (3o5), it is noticed, would permit reduction of the last three 

equations (4 0 1) to only two. However, the difference between the value of 

r 6 

jl obtained from the integrated values of v , v , v and that obtained directly 
by the Appleton-Hartree formula offers a method of checking the accuracy of 
the integrations at all points 0 

For the purpose of integration in the computer a seventh equation is 
added ^ 



(4.6) 


4.2 Initial Wave Normal Direction 

— — — — ■ “ r 9 '*<P 

If the initial values of r, 0, v , v , and v are known, the integration 

can be performed. In the case of a satellite transmitting and a ground 

station receiving, the values known are the position of the satellite 

(r 9 0 , *P ) which will be called point T, and the position of the station 
o o o 

(r , 0 , *P ) which will be called point R. 

111 r 0 <p 

A direct calculation of v . v , v at point T is difficult, but if the 

7 V. 

direction cosines a, B, Y of a vector in the initial ray direction referred 
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r 0 <P 

to the local orthogonal axis are known, v , v , and v can be calculated from 
them in the following way: 



(4.7) 


where v is related to the phase velocity in the direction of the ray by the 
PS 

velocity of light. 


... 1 dr 

v dt 

Pg 


= a 


1 d 6 

r — = S 

v dt 

Pg 


(4.8) 


1 . fi 

r sin o — = Y 

v dt 

Pg 


Substituting (4.4) into (4.7) 


v = 
Pg 


\j 1 + D 2 sin 2 © 
M- 


(4.9) 


r 0 C 0 

Substituting (4.9) and (4.4) into (4.8) and solving for v } v , and v : 


v r = H- 


a(l + D 2 sin 2 - D cos T 

i 

1 - D cos 0 


e 

v = h 1 


2 2 1 /2 

B(1 + D sin 0) - D cos "H c 


1 - D cos Q 


(4.10) 


<P 

v = H 


7(1 + D 2 sin 2 0) 1/2 - d cos n 


W 


1 - D cos 0 
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cos ©can be determined from (4 0 11) by an iterative method taking advantage 
of the fact that for frequencies fairly large cos ® = A 

After cos© and | sin ©I are known, v** 9 v and can be calculated from 

a,& and Y using Equations (4 o 10). 

4,3 Initial Ray Direction 

The problem now is how to calculate a , B and Y at point T in such a way 
th^t the ray will pass through point R. Since this solution is intended to 
be independent on the electron density distribution and the earth’s magnetic 
field models used, no mathematical derivation will be of help since it would 
be different for each mode 1 . In addition it would be extremely complex. 

A successive approximation method can be used. 

A ray integrated with the values a ± , at point T generally will 

reach the sphere r = r , with coordinated 9 and <P differing from those of 
point R, but at a point R^ * A correction applied to the values a ± , Y^ 

considering the result obtained, will yield the values a ^i+1* 
These values will be more successful in approaching point R 0 


This is explained in Figure 2a where: 

g s= unit vector in TR direction of components a , B , Y 
to o o o o 

g , g = unit vectors of components a., B>, Y. and a. ■, B. Y. 

& i > & i+l l* l i+l' i+l* l+l 

which would correspond to trials of order i and i+1 
g x = unit vector in the Tr[ direction with components a. . B^ ♦ Y^. 

&i i i* i’ i 

R / s the point at which the ray obtained with the values a., 8.. Y! 
i i' i' i 

reaches the condition r = r 

In this instance 9 a good correction when a plane case is considered, would 
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be to take g. as* shown in Figure 2, 

U \J 

g Q equal to the angle between g^ and g 


specifically, the angle between 

. ,, or if a veotor A is defined 
14*1' 


u/ 

g i 


and 




(4, 13) 


The correction explained above would make the points F and B symmetrical with 
points E and C in respect to the direction of the vector A. 

This latter criterion can be used in the three dimensional case to 
correct the initial ray direction used in trial of order i, as it is explained 
in Figure 2b„ The symmetry of F and E with respect to the direction of A 
can be expressed: 




= KA 


(4.14) 


1/ 

g i + l 



where K is a scalar constant of proportionality which can be calculated by 
substituting g ^ from the first of (4 0 «14) into the second: 

:• 2 gj . A 

K = — — - (4.15) 

A • A 

Substituting (4.13) into the first Equation (4,14) and expressing the equation 
■obtained, as well as (4*15), in terms of a ? s, B*s and Y 9 s, 


K = 


a\ ( a . + a ) 4- B / . (B. + B ) + V 7 OY. + V ) 

lio lio lio 

l + a. a +B, B + Y. Y 

1 O 1 o X o 


a. . = K (a. + a ) - a'. 

1+1 '1 O 1 


8 i + 1 ■= K < 8 i + 8 o> - 8 1 


(4.16) 


\ ♦ 1 - K <*i * V - * 1 
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The first trial can be the values a . B and \ . 

o' o o 

The following formulae, which give the components of the unit vector 
on the direction between two points T and R, at point T, can be easily 
derived from geometric consideration: 


r i\ r i 

d = / 1^+| — ] - 2 — [sin 6 sin 0 cos (<!P - ) + cos 0 cos 0 ] 

iT J r 1 o ol 1 o 

o' o 


a o d 


1 - — [sin 0_ sin 0 cos - <p ) + cos 0 cos 0 1 

r 1 o o 1 1 o 


(4.17) 


i = — \ • — [sin 0 cos 0 cos (<P - <P ) - cos 0 sin 0 ] 

o d r o o lo 1 o 

L o 


Y 


o 


1 

d 


r i 

[— sin 0 sin (<P 
r 1 1 

o 


V' 


Equations (4.17) are also used when calculating a / ., B / ., Y / . using the 

/ / ill 

appropiate values for 0 and . 

4.4 Coordinates 

The ability to give the coordinates of points along a ray, in order to 
plot a ray if desired, is a desirable feature of this program. This feature 
involves the capability of obtaining points along the ray referred to a convenient 
coordinate system. A local spherical system has been chosen, centered at the 
center of the Earth with the z / axis passing through the receiving station S 
and the x / z / plane passing through the satellite T. 

In this way r remains untransformed but the difference between r and the 
radius of the Earth a results in the height above the ground. The coordinate 
0 gives the angular distance from the station. The value gives the 
deviation of the ray from the vertical plane which contains T and R u Further- 
more, since the lateral deviation is small, multiplying ^ by a . sin 0 the 
deviation can be expressed in linear units perpendicular to the x / z / plane and 
measured on the surface of the Earth. 



20 


* It is convenient to use an intermediate coordinate system to transform 

t ) •*" // j / 

the coordinates,, This system is taken with the z axis coincided with 2 
and the x" z plane containing the x axis 0 In this case 


z 1 - z" 

Q / = Q" 
cp / = - cp 

o 


( 4 . 18 ) 


where is the coordinate p" of the satellite. The transformation from the 

o i. 

original to the intermediate coordinate system can be easily obtained. 


cos 0 = cos 6 cos 0^ 4* sin 0 sin 0^ cos \ 

1 1 1 \ 


tt sin 0 sin ) 

sin <P = — — 


in 0 


sin 


( 4 . 19 ) 


n cos 0 - cos 0 cos 0 

cos <p = 7/ 

sin 0 1 sin 0 
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5. DISCUSSION OF THE RAY TRACING EQUATIONS 


The differential equations used in the integration are the set (4.4), 

(4.6) and the last three Equations of (4.1). Interpreting the variables in 
these equations is somewhat confusing, particularly the variable t, therefore, 
a brief discussion, with a special case, may add understanding to its meaning. 

The point, A, and B shown in Figure 3 will be examined. Point A is 
associated with the set of values (r, 9, t) and point B is associated with 
the set of values (r + A r, 9 + A 9, + A % t + A t) . For discussion purposes 

let us choose point A on the ray such that H* is in the r-direction and the 
magnetic field is in the r 9 plane. These conditions are expressed as: 


* M- 


0 <p 

v = v = 0 

0 = T) 

r 

cos = 0 

Replacing the variable t by T: 


(5.1) 


t = c T (5.2) 

and taking finite differences in Equations (4.4) according to the conditions 
expressed in (5.1), the following results are obtained: 

* • * c 

Ar = jx A T 


‘ ■ r A 6 = 


c 

^2 80 


A T 


A <P = 0 


(5.3) 


The fact that point B belongs to the group trajectory can be checked 
by investigating the angle a between the wave normal and the obtained A s. 


tan a 


rA0 _ 1 
Ar “ H- 90 


(5.4) 

* 
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The last expression is a well-known relation between wave normal and 

10 . 

group velocity directions (Bremmer, 1949. ) . The first Equation (5.3) relates 
A r and AT by the familar phase velocity. The wave front AD will move to 
CB in the time AT, or the wave front moves a length As along the group 
trajectory in the time AT. This last interpretation and Equation (5.2) 
justifies the definition for the variable t. 

From Figure 3 it can £*lso be seen that the relations (2.5) are satisfied. 
The Faraday rotation cjin be expressed in terms of the variable t in 
the following way: / 

radians 

(5.5) 

turns 


to 

S2 = -T- (t - t ) 
2c o x 


— „ ( t — t ) 

2c o x 


The Snell law is the second point to be discussed. It has been shown 

that in a rectangular orthogonal system of reference, the derivatives — 

9|x ’ dx’ V 

-g- can be considered the components of a vector perpendicular to the stratification 

surfaces, and because of this property the rectangular coordinates will be 

used in this discussion, leaving to Appendix II the demonstration of the equivalence 

of the set of differential equations in both coordinate systems. 

The equivalent differential equations in rectangular coordinates to the 

last three Equations (4.1) are: 


dv X _ 3- 
dt “ H- 


dv y _ 1 fy- 
dt M- 3y 


dv Z _ 1 9|i 
dt ~ (-•- c$z 

Taking finite differences between two points it follows from (5.6) that 


x y z 

Av _ Av _ Av 

~ 3M- " SH- 

dx <?y ~&z 


(5.7) 
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Equation (5*7^ shows that A^ is parallel with the normal n^, to the stratification 
surfaces in the vicinity of the points considered* The situation is illustrated 
in Figure 4 * Snell’s law can be written^ 


Sin i l = ^2 Sin 1 2 


(5.8) 


It is seen that the stratification surfaces have been assumed constant 
between points 1 and 2, separated by As* The direction of the normal to 

a|x an 

this surface is determined by the value of the partial derivatives - 5 —, 

6|jl ^ . 8 y 

a nd evaluated at point 1, or between points 1 and 2. If a third point is 

considered, the partial derivatives will have, in general, a different value 

and the normal n between points 2 and 3 will be different from n . This is 
2 1 

why differential equations used in this method make no assumption as to the 
electron density distribution or the magnetic field model. In particular, 

no assumptions about horizontal or spherical stratification are needed. 

A 0.0 

A limitation in Equations 4„6 arises i when u approaches 0 , or 180 . 

d CD dv^ a 

— and — - tend to 00 , because of the presence of sin 0 in the denominator of 
dt dt 7 

these derivatives* The physical meaning is easy to understand because in a 
ray passing through the z, axis the quantity <P increases abruptly by IT* The 
derivative becomes infinite at these points* 
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Figure 4. Snell’s law between points 1 and 2 
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6 . GENERAL DESCRIPTION OF PROGRAM 


The general description of the program is given by the flow diagrams 

of Figures 5, 6, and 8 0 These diagrams will be explained in this section 

following a brief discussion concerning the integration method„ 

The differential ray equations are integrated using the ILLIAC Library 

11 12 

subroutines (Fl-Wheeler, 1953 ) and (F5-Muller, 1956 ), Both subroutines 

use the Runga-Kutta or Gill-Kutta method of integration of a system of 

13 14 

ordinary differential equations (Fosdick, 1958 and Gill, 1951 ) . They 

require an ’’auxiliary subroutine” which will calculate the following quantities: 


m dr d9 , m d^ 

2 dt At > 2 « At > 2 dt“ 


r 0 (p 

m dv m dv J m dv 

2 At > 2 777“ At, 2 At 

dt J dt 7 dt 


(6a) 


m dt m ^ 

2 it At = 2 At 


where 

At = Length of integration step 

2 m = A factor introduced to increase the accuracy of the integration. 

This factor is chosen considering the scaling of all the intervening quantities 
and the range of values that they can take. 

The difference between both subroutines is as follows: 

When an entry is made to the subroutine (F5) in a part of the program, 
the subroutine will integrate the equations until one of the variables reaches 
a given value before the control is returned to the main program. ° This means 
(F5) will integrate*- as many integral steps At as required. It will also 
determine the final fraction of At to adjust the given variable to the given 
value , 

Subroutine (F5) is used to integrate the diffential equations up to 
the surface of the Earth, which corresponds to a value of r equal to the 
radius of the Earth 0 
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When an entry is made to the subroutine (FI), the subroutine will 
integrate one integral step At and will return the control to the main program. 
This subroutine is used only when the coordinates of points along the ray are 
required. The last point does not concide in general with the point R. 

6.1 General Flow Diagram 

This diagram (Figure 5) shows the operation of input and the interpretation 
of orders in the data tape. 

When using the program, the first step is to read the program in the 
ILLIAC in the usual manner. The ILLIAC will stop in a "black switch stop" 

(B. S. S.). As soon as the desired portion of the data tape is set in the 
reader, the program will read the corresponding section. of data by moving the 
"Black switch up". Unless the data corresponding to the ray sought is the 
one. ini question, the program will stop again in a B. S. S. ready to read 
another section of data. The data of magnetic field, electron density, and 
frequency can be read in any order provided the proper order is heading 
the data. This will be explained in a future section. 

When the "ray sought data" is read, the program follows another pattern: 
it prints the coordinates of the points T and R and it begins to integrate the 
ordinary or extraordinary ray, as desired. After obtaining the extraordinary 
ray, the program is ready to read a new set of data. 

6.2 .Integration Flow Diagram 

Figure 6 shows how the program makes successive approximations of initial 

conditions. With the first trial values of the initial direction cosines 

a. 8, "V > the program calculates and integrates the corresponding values of 
$ 0° (p ° 

V v v'. After the integration is completed up to the surface of the Earth, 


<p -<p' 

and 0 - 9 9 

o o 

o o 


are small enough. 

If they are not it calculates a^, 6^ and repeats the cycle as many times 
as necessary. 

Figure 7 shows an example of the values W - ( P > I and 0 - 0 / versus 

-. .. 1 o i| o i| 

the trial number i. It can be seen that it follows an exponential law until 


- <P and 

n i 


e - e' 

o i. 


approach the limit of the computer's accuracy. 


In the exponential region we can write c 


e 

*5* 

CD 

1 


0 

CD 

I 

o 

1 


o 

O 

(p 

- c 0 / 


<p 

- <p 

o 

1 


o 



-v 


~B i. 

> 2 


( 6 . 2 ) 
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Figure 5. 
Note : 


DATA TAPE 

(EACH GROUP OF DATA 
IS PRECEDED BY 
AN ORDER) 


26 2I7N (order) 
FREQUENCY 

26 219 N (order) 


MAGNETIC FIELD 
PARAMETERS 


26 221 N (order) 


ELECTRON DENSITY 
PARAMETERS 


26 227N (order) 
RAY SOUGHT 
DATA 


(int)subroutine 


PERFORMS 

INTEGRATION 


INPUT SUBROUTINE 
(I N P) 


♦ READS FREQUENCY 
I 


READS MAGNETIC 
FIELD PARAMETERS 


♦ READS ELECTRON DEN- 
SITY PARAMETERS 

I 


♦ READS RAY SOUGHT 
DATA 

i 

PRINTS COORD. SATELLITE 
STATION 

IS RAY ORDINARY 
OR EXTRAORDINARY? 


MASTER ROUTINE 


SETS (INT) < 4 — 6- 


FOR ORDINARY RAY 




B.S.S. ( 24 185) 
SETS(INT) FOR«-d> 


B.S.S. 


(24 999 ) 


B.S.S. 


(24 999 ) 


B.S.S. 


(24 999 ) 


IF | IF 

ordinary extraordinary 


EXTRAORDINARY RAY 
I 


B.S.S. (24 181) 
SETS(INP) 


t 


NEXT ORDER FROM DATA TAPE 


General flow diagram 

B.S.S, means ’’black switch stop”: The ILLIAC stops. By moving 

up the ’’black Iswitch” it continues working. 
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Figure 6. Integration flow diagram 

Note: (P16) is a library subroutine for punching out numbers 














LOG (DEGREES) 


TRIAL NUMBER 
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where the constants B and B are proportional to the slope. of the 

-L ^ 

corresponding curves in Figure 8. The B’s depend on many factors, such as the 

distance between T and R, the electron density distribution, and to a lesser 

extent, the magnetic field variations between T and R. 

When the distance between R and R X is satisfactorily samll, the program 

will print the obtained value of t. If points along the ray are required 

r .0 

it will also integrate the ray again by using the initial values v ,v ,v 
used in the last trial. This integration will be carried out with the sub- 
routine (FI) and it will print out points along the ray and the corresponding 
values of the error. 


^ ■ i < - -V 


where H* = refractive index calculated from expression No. 3.5 using the 

« r 0 CP 

integrated values of v ; v } and v . 

\i 2 = refractive index calculated from Appleton-Hartree formula in the 

A .. : 0 cp 

direction defined by the values v , v , v used to calculate K 

It will be noted that subroutine (1-2) needs the use of the subroutines (ED) 

(EMF) and (NSQ) ; so does the subroutine (EWQ) . The last subroutine, is explained 

•in more detail in Figure 7. 

' 6.3 : Au*jl iary r Subroutine Flow-Diagram 

The function of this subroutine has been explained as complementary to 


(FI) and (F5) . 

The flow diagram corresponding to this subroutine is given in Figure 8. 

It can be seen that the subroutine (EWQ) requires the use of (ED) and (EMF) . 

(EWQ) calculates, then, cos ©and sin ©as the angle between the phase velocity 

r 0 (P 

direction defined by v , v and V and the magnetic field. With the values 

of cos © and sin © calculated in this way and those of N and H, (NSQ) calculates 

F- 2 and D using Equation Nos. (3.1) and (4.5). The final result of an entry to 

■ r 0 C Q 2 

subroutine (1-2) with the values of n, 0, Y, v , v and v is to obtain F 
o 

and D or F only from the local physical conditions. It should be noted 

r 0 (p 

that (1-2) uses the set v , v and v only as a direction. 

, 3F 1 3F . 1 9F . . . . 

The subroutine (DER) calculates gp , - |^j, and - -g ^ by taking 

finite differences. If the values AjT, A ^ and Zi <P are given: 
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Figure 8. Auxiliary subroutine flow diagram 
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fy; _ l_ :9(H- 2 ) 

9r 2H- §r 


u 2. V 

H- (r + — 


6 , <P) ~ 


H- 2 (r 


V 


2F 


9, % 


9jjl 

9B 


H- 2 (r^0 


V 


2r H- 


9 ) - 


lAi^S 




<P> 


(6.4) 


1 9H- 

r sin B cK p 


ix 2 (r e - 

v 1 1 ' 2 y 

'2H- r sin 


2 

JJ- (r Q 9 — ) 

^ v 1 l r 2 ' 




A^~ 


r 9 

where v , v , v remain constant,, 

This method of taking a derivative is equivalent to considering a second 
degree curve passing through every set of three points: (r + (A^iO/2, 0^ (f) 

(r, 0, <P) and (r - (A^r)/2), 0, <P) and similarly in the 0 and <jP directions 0 
The advantage of the method is that it can be used with any model of 
electron density and magnetic field a However, it has been found that to keep 
the error small, the second partial derivatives of N have to be continuous: 

This is depicted in Figure 9, where the error € (Equation No* (6*3) is shown 
versus height for a vertical ray in a spherically stratified isotropic 
medium with assumed "smooth" parabolic electron distribution* . The smooth 
parabola has been obtained by rounding off parabolic distribution with additional 
parabolas in the edges* This is done c . in such a way that the first derivative 
of N respect to h is continuous everywhere, but not the second derivative./ 

The error jumps at the points where the second derivative is discontinuous* 

In Figure 10 the variation of € is shown for a simple Chapman distribution 
of electron density* The corresponding errors are plotted versus height in 
Figures 10 and 11 for different values of A^r „ It can be seen that when 
AjT is relatively large, a systematic error is produced* Reducing the increment 
A-^r the systematic error decreases to a point where random errors appear. 

A further decrease of A^r results in an increase of the random error* The 
error € is not a function of /Vr alone; but it depends on many additional 




and « as a function of height 
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Figure 11. « as a function of height 
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factors which have been kept constant when the curves of Figures 8, 10, and 
11 were obtained. 

The subroutine (EWW) is a straightforward computation of the values (6.1) 
with the help of the ray tracing equations. 
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7. SUBROUTINES USED 


The theory and function of the principal subroutines have been explained 
in Sections 4 and 6. This section will be devoted to subroutines which have 
not been mentioned or which deserve special attention, specifically, the 
(NSQ), (ED), and (EMF) subroutines. 

7.1 ILLIAC Library Subroutines: 

15 

(Rl) Square root - Given x > 0 it will calculate ^x~ (Wheeler, 1958 ) 

X 10 

(S4) Exponential - Given -1 < x < 0 it will calculate e (Goldberg, 1956 ) 

17 

(T5) Sin-cosine - Given 0, it will calculate 1/2 sin Q (Wermer, 1958 ) 

(P16) Infraprint - Given x it will print correctly rounded up to 12 

figures in fractionary or integer form with the decimal point 

18 

anywhere (Gillies, 1958 ) 

(N12) Infraput - It will read sequence of numbers up to 12 figures in 


fractionary or integer form (Gillies, D. B., 1957 ) 

(Y5) Transfer block of words between the Williams Memory and the 

20 

Magnetic Drum (Elliott, 1960 ) 

Additional Subroutines: 

1 [l 2~ 

(SEC) Given x < - it calculates J — - x , or given 1/2 sin 0 it 


(RP) 


calculates 1/2 cos 0. 

o o 

Given 1/2 sin 0 it calculates - 90 < 0 < 90 


(SEC) and (RP) have standard entry, 

7.2 Special Subroutines: 

1) (NSQ) This program calculates the ordinary or extraordinary values 
2 2 

of H* and D or H- only if the following data has been conveniently prepared: 


0^ COS 0 

This subroutine has the following constants stored: 


H, N, sin W cos & and f. 


eH- 

K ° 

1 “ 2ff m 


K = 
o 


47T 2 m€ 


(7.1) 


(NSQ) calculates in the first seven words the parameters 



K 

oo 



K 

o 


(7.2) 
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and (NSQ) is arranged in such a way that, if the same frequency is used through- 
out a certain program (after the first entry to this subroutine) the next 
entries can be done to the right hand order of the word number seven of 
this subroutine thus saving computation time. 

A shortened version of the program has been prepared for the case in 
which the magnetic field is disregarded. In this case Equation 3.1 reduces 
to: 

l^ 2 = 1 - x (7.3) 

2) (EMF) Given r, 0 and this subroutine will calculate the earths 
magnetic field H and its components referred to the local orthogonal system. 
Three models have been used: 

a) Constant magnetic field 


H = const H Q = const 

r y 

b) Matched dipole approximation 



const 


H 

r 



[g° COS 0 + sin 0 (g 1 cos + hf sin <P) ] 
1 1 1 


(7.4) 


H fi = (— ) [-g° sin 0 + cos 6 (g 1 cos + h 1 sin <P) ] 


(7.5) 


[-gj sin ^ cos W 


The constants g^, g_^ and h* are to match the three components of the 
local magnetic field. 

c) Complete field model (Jones and Melotte, 1953 21 ) 
n+2 

H = - EE (— ) (n+1) H m (g m cos m ^ + h m sin m 

r n m r n & n n T/ 


n+2 

=-cosec 0 EE (— ) m H m (-g™ sin m <P + h m cos m <P) 

0 nmr n & n n 


(7.6) 
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=-cosec 9 


n+2 dH m 

E Li (—) sin 0 a 
n m r do 


(g m cos m <P + h 1 " sin m <P) 
n n 


m 


(7.6) 


Where : 


H (cos 0) 
n 


2^ n ! (n - m) ! 
(2 - n) ! 


P (cos 0) 
mn 


(7.7) 


P = associate Legendre polynomial 

mn 

g , h .= coefficients with the dimension of amper per meter. Those 

n n ™ 
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calculated by Finch and Leaton (1957 ) are tabulated in 

Appendix I, 

3) (ED) This program calculates the electron density as a function 
of the coordinates (r, 0, <P) 

a) Simple parabolic layer 




N = 0 


for r - r I < A 


for r - r > A 
m 


(7.8) 


where 

N = maximum electron density of the layer 
m 

r = value of r corresponding to maximum electron density 
m 

A = semi-thickness of the layer 
b) "Smooth" parabolic layer 


N = N 


m 


1 - 


(r - r ) 
m 

ca 


2 1 


for 


r - r 


ml 


< c 


N = N 


(a - r - r )' 
m 

m a(a - c) 


for c < 


r - r 


ml 


< a (7.9) 


N b 0 


for 


r 


r > a 

m — 


where 

a = semi-thickness of "smooth" parabolic layer 

A = /ac = semi-thickness of equivalent simple parabolic layer 
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c) Simple Chapman layer 


-(r - 


N = N e 
m 


H - (r-r) -He 
m 

2H 


H 


where 


H = scale height of the layer 

d) Chapman layer with a gradient along meridian 


-<r - r ) 
m 


N = N e 
m 


H- (r-r) -He 
m 

2H 


H 


(7.10) 


N (.0 - 0.) + N . (0. 

_ m, l + 1 ■ i m, i » l + l 

m “ 0 ~ - 0. 

l+l l 


" (Q - 6. ) + r . . (0. 

m, i+l’ l m, l i+l 


0) 


- 0 ) 


(7.11) 


r = 
m 


“0 r-r 

i + l i 


H = 


H 1 + 1 • <e - V * H i • <9 i + 1 - e > 


■B r-g- 

i+l i 


where 

9 i + 1 > 6 i for i = 0, 1, 2, 3, A, 5, 6, 7 

The values N r H. corresponding to the values 6. are stored in 
mi' mi 1 1 


the program. 
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8. SCALING AND LIMITATIONS 


The ray tracing differential equations have the following limitations: 

a) Limitations inherent in the ray theory discussed in Section 2. 

b) Limitations in the value of 0 mentioned in Section 5. 

c) Limitations to be explained in this section concerning the range of 
values that the intervening variables are allowed to take. 

These additional limitations ^ arise from the following two facts: 

(1) The computer will not operate with numbers equal to or larger than one. 

(2) The necessary scaling of the different variables in a convenient way 
so that they will accommodate the requirements of the computer. 

The scaling of the variables is difficult because they have a very 
large range of variation,, Another difficulty arises because accuracy is 
lost when full use is not made of the available digits in the computer. 

A compromise has been made in this particular program that keeps in mind 
its application to the use of data collected from artificial satellites. 

A floating point arithmetic has not been used because it would require 
exceedingly long computer time. This makes its use prohibitive. 

The scaling of the principal variables is shown in the following list; 
where the primed variables are the scaled ones: 


a 


d / 


= r (Kin) X 10 

= t (Km) X 10 

_ a (radians) 
77 

a (degrees) 
180 

= d (Km) X 10 


a is any angle 
d is any distance : 


( 8 . 1 ) 


v 


X 2 


-2 


' -2 
^ = M* X 2 

_ JL 

// 2 
^ = F* X 2 


i = r, Q, (p 

two scalings are used 
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sin a = sin a X 2 
D* = D x 2 _1 


a is any angle 


/ -14 

N = N (electrons per cubic meter) X 10 


k f —2 

= H (ampers per meter) x 10 

/ -39 

f = f (cycles per second) x 2 


In addition to the direct scaling of the variables, the scaling of simple 
functions should be considered in the Auxiliar Subroutine and in (NSQ) . 
The following simple functions are defined in the auxiliar subroutine: 


1 — x [v r - 4D / (v r cos' 0 - p. / cos * *1 )] 

/ a T* 


1 r 9 / , 0 f / 

— j ~ 2 l v “ 40 ( v cos ^ - H* cos Tq) ] 


1 r y' ,y 


[v r - 4D (v r cos cos' t| ■)] 


in 7 0 


Ayl 1 2 = H- //2 (r + ^ AjT, 6, <P, v r , v 6 , v^) - H* 2 (r - i A^r., 0, <P, v r , v^) 


,/2 = H- /2 (r, 0, + i A x 0, ( P, v r , v^) - p.' /2 (r, 0 - A^, 9, v* , V Q , 


, r „ 0 „"2 « 1 fl r .0 


a/' 2 = r" 2 (r, 0, <P.\ ttf, v r , v G , v*> - (or, 0, <P - \ A/, /, V G , v*) 
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/ //2 

A r = 8 A x r V- 


An = 


A,„ = 




// 2 / , 
8tt H* r A i 0 / 


16tt r 7 sin 7 0 


( 8 . 2 ) 


At = 2 m At' 
m 


The values (6.1) obtained from the ray tracing equations can be written with 
the new scaled variables in the following way: 


m / dr 

2At dF 


2E A- t 
r m 


m / d0 
2 At 

dt' 


/ 

r t r 


V 't 
m 


m / d^ 

2 At 

/ 

dt 


A t 

/ / n Y m 

r ir sin . ti 


(8.3) 


2 m A t' HL 
dt' 


A/ 

= <-A— + K e E 0 + WV 

r 


m / / dv0 

2 A t — : — 


dt 


V' 


- ( -ir- - V, * 3 V V V 


//2 

SV' — = - y r " B K ( ^ e )A m t 


dt 




Each of the functions defined in (8.2) is independently scaled in different 
parts of the program. Each of them also has its own limitation of not 
exceeding the value one. 


E = E 
r r 


E e = E e 


(8.4) 


® cp “ E cp 
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K Q = K 0 X 2 


-4 


K '(p - K^, X 2 


/ —4 

B = B x 2 


a' i* 2 = a V-" 2 X 2 7 

r r 


/ 2 //2 7 

A qH- = AqH- X 2 


/ 2 // 2 7 

“ V X 2 


A' = A x 2 16 

r r 


/ 16 
A Q = Ag X 2 


/ . _16 
A <£> “ A <£> x 2 


A 7 t = A t x 2" 

m m 


(8.4) 


The scaling of the simple functions defined in (7.2) which occur in 
(NSQ) are: 


k' = K 2 -3 x 10 14 
oo oo 


K 'n = K n x lo2 


(8.5) 


The range of values that the different variables and other intervening 
quantitites can take in an application of the program can be deduced from the 
scaling (8 . 1) , (8.4) and (8.5). In each case, the most severe limitations 
should be considered. 
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The only exception which can be listed, of a variable exceeding the 
value one without affecting the end result, is that of scaled angles. When 
the value of a scaled angle exceeds the value one, it takes the value of the 
equivalent angle with the sign changed. The following example of the ILLIAC 
arithmetic explains this feature: 

175° + 20° = - 165° (8.6) 

The range of allowed values of each variable will be examined in detail. 

Electron Density - From the tenth equation of (8.1) 

N < 10 14 (8.7) 

Magnetic Field - From the eleventh equation of . (8 «, 1) 

H < 100 Ampers per meter (8.8) 

The scaling of the electron density and magnetic field satisfies the 
extreme conditions which occur in the ionosphere and the earth’s magnetic 
field. 

Frequency - From (8.5) 

3.52 me < f < 254 me (8.9) 

Refractive Index - From the seven scaling of (8.1) and the first three 
of (8.4) 

-1 < H - 2 < 2 (8.10) 

4 

The maximum limit is 2, however, the region of interest is: 

I < H- 2 < 1 (8.11) 

When Y is small, (8.11) corresponds approximately to: 

| > X > 0 (8.12) 

4 — 

r 0 (p 

The values v , v and v are components of a vector of magnitude and the 
only limitation is: 

v* < 2 for i = r, 0 , <P 
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Finite Differences A^r, £ Q and A 

These values are constants stored in the subroutine ( DER) , and could be 
changed if desired, providing that the limitations (8.1) are considered; 
otherwise, the corresponding scaling should be changed. From the seventh, 
eighth, and ninth of (8.4): 

< • 3814 Km 


A i« < 7 


3814 radians 


(8.14) 


A^ < 


1 

r 


3814 radians 


The stored values are: 


A^r = .22 Km 

& Q = 1.5708 X 10~ 5 radians (8.15) 


= 1.5708 X 10 


f adians 


Partial Derivatives of Refractive Index 

Two limitations apply to the derivatives: The first one concerns the 

finite differences of the refractive index and their own scaling given in 


(8.4): 


A M* < 2 


-6 


2 -6 
A e H- < 2 


(8.16) 


v ** 2 < 2 ~ 6 

The other, and more significative, limitation arises from the ratio of 
the scaled values which appear in the first term of the right hand of the 
last three Equations (8.3): 
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^<f ! x* 1! x 10-V 1 

Or i ' 


9(u 2 ) 2 12 -5 

gg - < r V- X 2 X 10 


(8.17) 


^< IS i.9 f 2 2 12 X ID’ 5 


Distances 

In general, the distances have the only limitation which arises from 
the corresponding scaling: 


d < 100,000 Km (8.18) 

Coordinates 

The limitations for r are from (8.14), (8,15) and the fifth (8.4): 

5,250 Km < r < 24,280 Km (8.19) 

r could take values larger than the one mentioned if the values (8.15) were 
chosen differently. 

The limitations for 0 come from the scaling of S given in the sixth . 
of (8.4) 

3° 34* 35" < 9 < 176° 25 7 25 V (8.20) 

No limitation exists for <P due to the arithmetic of the ILLIAC for the 
angles shown in (8.6) 

Variable t - Its limitation arises from the scaling given in (8.1) 

t < 100,000 Km (8.21) 

Integral Step Length - From the last of (8.4) 

m ■ q ^ 

2 At < 2 10 = 160 Km (8.22) 

Since m > 1, the maximum integral step length is: 

At = 80 Km 


f 

// 


(8.23) 
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9. PREPARATION OF DATA TAPE AND USE OF PROGRAM 


A data tape is prepared for each particular program. The flow diagram 
of Figure 5 shows clearly the different part of the data tape and the need 
for a heading order for each group of data. Each group of data will be 
examined in detail in the first part of this section. 

9.1 Frequency Data 

26 217 N Heading Order 

f Frequency in cycles in integer form 

J Terminating symbol which can be either of the following 

letters: N, J , F or L, 

9.2 Electron Density Data 

These data are prepared differently for different electron models. A 
simple spherically stratified Chapman layer is the only model used in the 
Faraday rotation program. 


26 221 N 

/ 

H 

r / 
m 


Heading order 


-5, 


N 


m 


Scaled scale height in (Km X 10 ) 

Scaled radius corresponding to the point of maximum electron 

✓ -5, 

density (Km X 10 ) 

Scaled maximum electron density of the layer (electrons 
per m X 10 ) 


J Terminating symbol which can be either N, J, F or L. 

9.3 Earth Magnetic Field Model 

These data are prepared according to the model used. , At present, two 

models have been used: 

Complete Field Model 

26 219 N Heading order 
1 ✓ 


1 

i/ 


49 scaled coefficients (amper per meter X Jqq') 


i/ 
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2 

h 




g 


o / 

6 


6 / 

V 

j Terminating symbol, either N, J, F or L 

Matched Dipole 

Heading * order 

0 , 2 V 

Scaled g (amper per meter X ^qq ' 

1 4 
Scaled g (amper per meter X ) 

i 4 

Scaled h (amper per meter X *Joo^ 

Terminating symbol, either N, J, for L 

9.4. Ray Sought Data 

This datum is formed by two parts: the fractionary numbers and the 
integer numbers: 

26 227 N Heading order 


26 219 N 
o / 


1 

1 / 
1 

1 / 
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Fractionary Data 




^ Scaled coordinates of point T or transmitter. 


/ 


<P, 


y Scaled coordinates of point R or receiver. 


m ✓ 
2 At 


Scaled value of integral step, 


Terminating symbol, either N,;J, F or L 


Integer Data 


m 


2 m is the scaling factor used in subroutine (F5) and (FI) 


a 



b 


c 


J 


If a = o, the program will use \he subroutine (F5) only.' 

If it is a positive number larger than one, the program 
will use the subroutine (FI) and will print the coordinates 
of points every ”a" integral steps. 

12 1 

The program will print ‘-H* and — € every b printed points, 

^ A . 12 1 

or every a.b integral steps. If b = 0, — and — € will 
not be printed at all. . ; 

If c = 0 the program will first integrate the ordinary ray. 
If c = 1 the ray will integrate the extraordinary mode 
directly. 

Terminating symbol, either N, J, F or L 


52 


All the numbers of the data, fractionary or integers, can be written to 
twelve figures . These numbers are preceded by the sign + or - . (Gillies, 
1957) 

The program has in itself values of data corresponding to frequency, 
electron density, and magnetic field models. The arbitrarily chosen written- 
in values ares 

For frequency: f =.20 me 

For electron density: H = 60 km 

(Simple Chapman layer* only): r = 6721 km 

m 

N = O 0079 

For magnetic field: The coefficients are those given by Finch and 

Leaton (1957) when the complete model is used. The coefficients for the 
matched dipole model are those obtained by matching the field with the 
"complete model" at 300 km above the receiving station at the University of 

Illinois. : 

Use of the program is fairly well explained in Figure 5. Listed below 

is a summary of the steps to follow: 

Step 1. Read the program in the usual fashion. I should stop in a 

24 999 (24 3F7 ) order at the right hand order of location 87 (057 ) . 

16 lo 

Step 2. Set the data of frequency, magnetic field, and electron density 

in reader* This can follow, any order. Move black switch once for each 

group of data. After reading each group of data, the program will stop in 

a 24 999(24 3F7 „) order at locations 218 (0 JK ) , 220 (0JN ) and 222(0JF ) 

16 16 7 16 lb 

respectively. If it is desirable to use written-in values, the corresponding 
group^ of data should be omitted. 

Step 3c Set the ray sought data in the reader and move the black switch 

up once again. This will cause the program to start working. 

Step 4c If the data has been prepared to integrate the ordinary mode, 

the program will stop in the B.S.S. 24 185 (24059^) at the right hand order 

of location 180 (0S4 „) when the ray is completed. Move the black switch 

16 

up if the extraordinary ray is desired. Return to Step 1 if a different ray 
is desired* 

Step 5. When the extraordinary ray has been completed, the program 

will stop with the B.S.S. order 24 180 (24 0S5) at the left hand of 

location 180 (0S4 „) . The program is then ready to proceed with a new ray. 

16 

Return to Step 2. 
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10. PRELIMINARY RESULTS 

The analysis of a passage of the satellite 1958A^ recorded at the receiving 
station of the University of Illinois was undertaken for the purpose of studying 
the behavior of the program described in this report. 

The passage recorded between 0900 and 0908 C.S.T. on January 11, 1960 
was chosen in accordance with the following: 

a) It was a passage close to the receiving station. The point 
of closest approach corresponds to a ground distance of about 
200 km. 

b) It was a low passage. The satellite was at or below the height 
of maximum electron density of the ionosphere. 

c) There was an ionogram sounding taken at the University of 
Illinois at 0900 C.S.T. 

The coordinates of the satellite were taken from the orbital data provided 

by the Smithsonian Astrophysical Observatory. These coordinates are shown in 

Figure 12. The points used for analysis correspond to times of integer minute 

and the points corresponding to time 0903, 0904, 0905, and 0906 will be called. 

Points 1, 2 , 3, and 4 respectively. Figure 12 also shows the M.U.F. predictions 

23 

published, by the National Bureau of Standards (1959) corresponding to the 
month of January. 

: k. The electron density distribution used was a spherically stratified Chapman 

layer (Equation 7 . 10) . This layer is an approximation to the profile deduced 

from the mentioned ionogram. A comparison between both profiles is shown 

in Figure 13. The ionogram was reduced by the method described by Wright 
24 

and Norton (1959) . The obtained parameters of the Chapman layer are: 

H = 62.7 km 

h = 300 km 
M 

12 3 

N = 1,3416 X 10 elec per m 
M 

The magnetic field used was a centered dipole (Equation 7.5) matched to the 
"complete field model" (Equation 7.6). This complete field mode uses the 
coefficients given by Finch and Leaton (1957) at a point 300 km above the 







HEIGHT Km 


ELECTRON DENSITY (ELECTRONS PER m 3 ) 

Figure 13. Assumed Chapman layer and measured electron density distribution as a function 
of height 



Ui 
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receiving station. The values of the magnetic field at the matching point 
are: 

= -37.5208 Amp. per m 
Hq = -12.9927 Amp. per m 
H ^ = + 1.0483 Amp. per m 
| H | = 39 . 7205 Amp. per m 

The ordinary mode rays obtained are shown in Figure 14. These rays have been 
plotted to gain an idea of the magnitude of the blending of the rays. 

The lateral ground deviation of ordinary and extraordinary rays are shown 
in Figure 15. 

m The Faraday rotation was calculated for the four points with the values of 

integral step At = 80 km. At = 40 km and for Point 3 with At = 5 km. Table 

1-shows the values of t,. t, x and £2 obtained in these calculations. It 

(oV (x) 

is interesting to notice that the resultant values of £2 are almost independent 

of; the values of At, since the variations of t, N and t, x are in the same 

’ (o) (x) 

direction and compensate each other. This observation should be taken with 
reservations when very low rays are considered since in this condition a very 
•' small error in the components of H- around the height of maximum electron density 
results in large difference on the path of the ray. 

The r.m.s. values of €, defined by Equation 6.3, have been calculated 
for the. ordinary ray corresponding to these points. Together with | € | 
these values are plotted as a function of phase length (t^) in Figure 16. 

Contrary to expectations, the shortest ray;- corresponding to Point 3, has 
the largest m.r.s. error. The same values of € are plotted as a function of 
azimuthal angle, measured clockwise from the North, in Figure 17. The 
resulting curves of Figure 16 and 17 are smooth even though the sequence of 
points is different. This seems reasonable since the variation of the error 
€ depends on other factors such as electron density distribution, (as will be 
shown later), slope of the ray, and the angle between* the direction of propagation 
and the magnetic field. It would be necessary to have a large number of 
systematic calculations to individualize the effect of each variable. 



height (Km) 



LATERAL DEVIATION IN Km 







TABLE I 


Coordinates of satellite 1958^ on January IT y 1960 and computed values of t 
t(x) and ^ with different integral steps- 


(o) j 


Point 1 


Point 2 


Point 3 


Point 4 


Time (C. S. T.) 
height (km) 
latitude (degrees) 

j 

longitude (degrees) 


09 hrs. 03 min. 
1,291.7 


46.88°N 


91 ,97°W 


09 hrs. 04 min. 
276.6 


43.63°N 


88.62°W 


09 hrs. 05 min. 
262.3 


40.27°N 


85 . 79°W 


09 hrs. 06 min. 
249.0 


36.82°N 


83 ,24°W 


At = 80 km 


At = 20 -tan 


At = 5 km 


t (o) (km) 

839 , 

.617 

440 

2 

479, 

.752 

390 

3 

334, 

,788 

727 

9 

617, 

.033 

499 

t (x) (km ) 

838 . 

.470 

507 

7 

479, 

.187 

669 

6 

333, 

.896 

402 

1 

615, 

.824 

195 

(rotations) 

38, 

.564 

42 


18, 

.824 

02 


29, 

.744 

19 


" 40, 

.310 

16 

W * 0 

839, 

.613 

209 

8 

479, 

.747 

446 

8 

334, 

.783 

677 

0 

617, 

.029 

795 

t (km) 

(x) 

838 , 

.462 

668 

,4 

479, 

.178 

110 

8 

333, 

.889 

749 

5 

615, 

.819 

564 

Q, (rotations) 

38, 

.351 

35 


18, 

.977 

9 


29, 

.797 

58 


40, 

.341 

06 

t (o) (lan) 









334, 

.783 

018 

3 











- 


333, 

.889 

701 

1 




S2 (rotations) 









29, 

.777 

24 






oi 

to 
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The simplest example of individualizing the effect of one variable, and 
a very interesting one in the behavior of the program, . is the following: The 

ordinary ray corresponding to Point 3 has been integrated using the values of 
At equal to 80, 40, 20, 10 and 5 km* The r.m.s. values of £ are plotted as 
a function of At in Figure 18. Note that the upper part of the curve between 
At = 80 km and At = 20 km follows the law: 

• ' e = a (At) b (10.1) 

r.m.s. 

where a and b are constants which can be calculated from Figure 18. .The' 
lower part of the curve decreases more slowly and it is expected that it will 
reverse its slope for smaller values of At, following a pattern very much like 
that explained in part 6, when discussing the variation of A^r . 

Decreasing the integral step At by the factor 1/2 increases the computation 
time by the factor 2. A value of At = 20 km is considered by the author a good 
compromise between computation time and obtained error. However, for rays not 
very critical, At could be increased to 80 km without appreciable difference 
in the final value of Faraday rotation. This has been pointed out with the 
aid of Table 1 . 

Another interesting feature is" noticed by plotting JE as a function of 

height (t, x or t, x could also be used to this purpose) as it was done in 
(o) (x) 1 

Figure 19 for the Point 3 and for different values of At. Variations of € 
occur only where there is appreciable electron density and remain constant 
below 80 km, where the electron density is negligible. This occurrence is 
more likely to be due to the presence of electron gradient rather than to the 
presence of the electron density alone. 

Figure 20 is a comparison of computed and measured values of the "Faraday 
rotation", from. the corresponding record obtained at the receiving station of 
the University of Illinois. The origin of the "Faraday rotation" axis for the 
measured curve \s arbitrary. It was brought to coincide to the calculated one 
at the point of closest approach. . The difference between ;the two curves is very 
large. However, an analysis of the assumptions used in the electron density and 
magnetic field models show that the difference is reasonable. The predictions 
of M.U.F, shown in Figure 12 indicate that the passage was at a time of the day 


o o o o o 

in O c\j ro 10 cc 



Figure 18. Variation of * 

& r.m.s. 


as. a function of At for point III 
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when the horizontal gradient is very large. A smaller electron density at the 
N.W„ would result in a smaller number of .rotations at the beginning of the 
passage. A larger electron density at the S„E. would result in a larger number 
of rotations at the last part of the passage. In addition ot this, it can 
be taken into account how critical the position of the satellite is to the 
electron density model used. The satellite was slightly below the peak of 
the ionosphere, and a small horizontal variation in the values of H and h^, 

added to the variation of N could very deeply affect the final results. 

M 

The northern part of the passage is also very sensitive to the magnetic 
field model used as noted by Little, hnd Lawrence, (I960). ' . ' . 

From the foregoing discussion it can be seen how useful it would be to 
have an improved electron density model to arrive at a method of calculating 
the electron distribution from the Faraday rotation using artificial satellites. 
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11. PROPAGATION PROGRAM 

A modified program, which considers a ray propagating in a meridian plane, 
has been prepared in order to investigate long range propagation problems. 

This program has the following characteristics: 

1. The ray is contained in a meridian plane. The coordinate system 
used is a spherical one centered at the center of the Earth. Its orientation 
is determined for each problem so that the ray of interest is in a meridian 
plane. 

2. The magnetic field is disregarded. Therefore, the refractive index 
is calculated with Expression (7.3). 

3. The electron density distribution used is defined by Expression (7.11). 

4. A downward going ray is mirror reflected when it reaches the surface 
of the Earth. 

5. The output tape is suitable for use in the data plotter at the 

25 

Digital Computer Laboratory of the University of Illinois (Carter, 1958 ). 

6. The allowed ranges of variation for r and 0 are given by (8.19) 
and (8 . 20) . 

7. The integration of the ray is interrupted either when r exceeds 

a given value r p or when 6 exceeds a given value 0 p , whichever happens first. 
Data Tape This tape consists of two groups of data: 

Electron Density Data This group of data consists of thirty-two 
decimal numbers up to twelve figures preceded by the sign plus: 

- Q* eight scaled angular distances 

- 0 corresponding to points for which ionosphere data is 

o 

available. 6 is the inivalue of 0 for the ray sought. 

/ 

- 0 

o 
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✓ 

H 


o 


H 


1 


scaled height at the points 
defined by the previous set 


H, 

N 


✓ 


✓ 



scaled maximum electron density at the same points 




scaled value of r at which the electron density is maximum 
at the same points„ 



N Terminating symbol, either N, T, F or L 

Ray Sought Data Consists of a set of fractionary numbers and a set 
of integer numbers. Each number is written up to twelve figures and preceded 
by a sign. 

24 99 N Heading order 

Fractionary set of data 

i Scaled initial inclination angle 0 This angle determines 

o 

the initial inclination of the ray, and it is measured 
from the up going verticle at the initial point e 
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2 m At 7 


F5 


F5 


N 

Integer set of data 
a 


b 


m 

N 


Scaled initial 0 
Scaled initial r 
Scaled final 0 
Scaled final r 

Scaled integral step. It must satisfy (8.22) 

Scaled r for (F5) test. 

Scaled v 1 * for (F5) test. r„ _ and. v** should be 
F5 F5 F5 

cho sen so that a down coming ray that reaches the 

r r 

height corresponding to r__ with a v smaller. than v __ 

F5 F5 

will hit the surface of the Earth. The effect of the 
(F5) test is to transfer the integration to the subroutine 
(F5) which will integrate the ray exactly up to the 
surface of the Earth without any output in between. The 

ray is reflected at this point by reversing the sign of 

r r 

v . Suggested values for r and v are .643471 and 

F5 F5 

- .0351 respectively. 

Terminating symbol. Either N, T, F, or L. 

Output frequency. The program will print out the 

coordinates of a point every a integral steps. 

Error output frequency. The program will print out 
2 

1/2 and 1/2 6 every a.b integral steps. 

2 m is the scaling factor used by subroutines (FI) and 
(F5). 

Terminating symbol. Either N, T, F or L. 


Use of Program The following steps should be observed when using the program. 

Step 1. Read the program in the usual fashion. It will stop in a 24 999 

(24 3F7 ) . 

16 

Step 2. Move the black switch up once. This will cause the ILLIAC to read 

an "Auxiliar program" at the end of the tape. The program will stop in a 

24 10 (24 K 1C ). 

16 
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Step 3. Put the electron Density Data Tape in the reader and move the 

black switch up. When these data are read, the program will stop in a 24 99 

(24 3F7,„) order at location 11 (S.,„) . 

16 16 

Step 4 0 Put the ray sought data tape in the reader and move the black 

switch up once. The program will start working* It will stop in a 34 999 

(34 3F7 ) order at location 126(07F‘* ) or 128(080 ) when the ray is completed. 

16 16 16 

At this point it is possible to go to step 2, or, if the same electron density 

data is to be used for a new ray, go back to Step 4. 

Example. An application of this program is shown in Figure 21. The propagation 
of two rays of 15 me between Slough, United Kingdom and San Francisco, California 
have been obtained using ionosphere data furnished by I.G 0 Y., U. S. World 
Data Center. The data are from stations near the meridian plane joining the 
two mentioned points and it corresponds to August 23, 1956, 00 hrs . 00 min. 
Universal time. 

The largest sources of error in this application originate from the 
difficulty in analyzing the available ionograms and also in assuming that 
the stations from which ionosphere data was obtained were located on the 
meridian plane mentioned. 



0001 
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APPENDIX I 


The constants used in the Faraday Rotation program are the following 
electron charge: 

“19 26 

e = 1.6020 X 10 coulombs (Cohen, 1955 ) 

electron mass: 

-13 27 

m = 9.10721 X 10 kg (DuMond and Cohen, 1951 ) 

velocity of light: 

9 -1 28 

c = <,2997925 X 10 m sec (Froome, 1958 ) 

magnetic permeability constant for free space: 

-7 

= 4tt X 10 henry per meter 
o 

dielectric constant for free space: 

£ = ■ farads per meter 

° i/ 2 

V 

Special constants for (NSQ) subroutine; 
eH 

K = —2 = 35180.9170975 M.K.S.(r) 

1 2irm 

2 

K = — ? — - = 80.6178714557 M.K.S.(r) 

2 4irra€ 

o 

Mean radius of the Earth: 

29 

a = 6371 km (Sitter, 1938 ) 

Mathematical constants: 

Tt = 3 o 14159265358979324 

e = 2 . 718281828459045 



75 


Coefficients for "complete 11 earth magnetic field (Finch and Leaton, 1957) 
are given in Table II 

Location of receiving station of the University of Illinois: 
lat: 40° 01 T 04.6" N 


long: 


88° 19’ 37.2" W 



TABLE II 


Gaussian coefficients in e.m.u system 



0 


+ •30545 
+ .02280 
- .02950 
-.04170 
+.02035 
-.01495 


+.02265 

-.05245 

+.05855 

-.04400 

-.03235 

-.01005 



-.05910 

+.03295 

+.01390 

-.00810 

-.00185 

+.00495 


2 

3 

4 

5 

6 

-.01370 





-.02445 

-.00720 




-.02260 

+.00795 

-.00230 



-.01550 

+.00190 

+.0034 

+.00050 


-.00190 

-.02350 

+.00165 

-.00010 

+.00070 

h m 

n 

2 

3 

4 

5 

6 

-.00215 





-.00565 

+.00080 




+.01185 

+.00080 

+.00120 



-.00715 

+.00210 

+.00315 

-.00060 

- 

-.01715 

.00000 

+.00055 

+.00060 

+.00010 



OJ 



































77 


APPENDIX II 


The purpose of this appendix is to sketch the transformation of Haselgrove f s 
ray tracing equations from the rectangular coordinate system of reference R to 
a spherical system, S. The method will be a straightforward I transformation of 
variables. The direct method has been chosen to avbid the use of tensorial 
calculus for the benefit of those who are not familiar with this branch of 
mathematics. In addition the tra;nsf ormation helps to increase familiarity 
with the meaning' of the various derivatives and also proves how much work is 
saved when tensorial methods are used in the transformation of these equations. 

The Haselgrove^ ray tracing equations in a rectangular coordinate 
system are: . / * 


dx 

dt 



A|jl . 
H* — ) 

3v X 



M- 


*-> 

d v y 




dv. x _ 1 ^ 3jx 
dt H- ?x 


dv y _ 1 9^ 
dt _ M- "3y 


dv^ _ ^ 3p- 
dt “ '9z 


* 


UI-1) 


From Figure 22, the variables of the S system can be written in function 
of the R system variables: 


r 


f 2 2 2 

Yx + y + z 


( I 1-2) 





x,y,z,r,0,<#> ARE UNIT VECTORS 


Figure 22, Relation between "R" and "s" variables 
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0 = tang 


-1 


J77 


<P = tang ^ — 


r 1 

v = 


. f~2 2 2 

Vx + y + z 


, . x y z \ 

(v x + v y + v z) 


0 

v = 


.f~2 2 2 f~2 2 * 

yx + y + z Vx + y 


r x y z 2 2., 

[v xz + v yz - v (x + y )] 




fFl 


, x y x 
C-v y + v x) 


2 

y 


The inverse transformation will be: 
x = r sin 0 cos <P 

y = r sin 0 sin 

\ 

z = r cos 0 

v x = v r sin 6 cos <P + cos 6 cos <P - v^ sin 

y r 0 

v y = v sin 0 sin + v cos 0 sin ^ + v cos 

z r n 0 . q 

v = v cos 0 - v sin ti 

Taking the total derivatives of (II-2) with respect to t, the result 
will be the sought derivatives of the variables of the S system in terms 
of the derivative of the R system in respect to time. Substituting the 

derivatives of the R variables by the values (II-l) and the R variables by 

\ 

the values (II-3), the resulting equations are: 


( II-2) 


(II-3) 
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•^7 = - i [sin 9 cos <P — — + sin & sin <P + cos 0 — 

dt - 2 11 9 V X a v y 9 v z 


a v y 


dO v 6 1 , fi ^ ft ■ co ^ ■ ft ^ l 

— = — zr - — [cos 0 cos Sr + cos 0 sin Y - sin 0 1 

dt ..2 rH- a x - ” - - J 

rH- °v 


a v y 0v 5 


I? ^-5 ( lf-i[ - , in <P + CO. 

dt r sin 0 |^2 H dy x dy y J 


dv r _ 1 , 0 2 <?\ 1 / , 9 0 m <P . rCK ^ 

— - (v + v ) < ( v cos y COS Y - V sin Y) + 

dt r|i 2 r| ^ l 0v x 


0 (p 3jjl 0 9[i I 

(v cos 0 sin <P + v cos <P) - v sin 0 r + 


( II- 4 ) 


Sv y 


■'V 


1 

+ jX 


d\i 9jjl 9jjl 1 

sin 0 cos -g^ + sin 0 sin <P -g^ + cos 0 -g^ r 


dv 1 


dt 2 

rH* 


0 $ cos 0. 


- V v + v — n) + 77“ i 

sin t) H-r 


A / r tn *P sin ^P\ 

cos o (v cos Y + v q ) 


sin 


q , r o m <P cos <p^ r Q fy- 

a v x sin 9 a v y 


9 v Z 


} 


+ i cos 0 cos ^ -g^ + cos 0 sin <P -g^ - sin 0 j* 
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dv^ 1 , r <P . a 9 <P Qx 1 f / r . q 

_ = (v v sin 6 + v v cos 0) + pr ri r n ~ e \ ~ (v sin 6 

r sin 6 l 


0 0|i 0 

sin <P + v cos 0 sin <P) + (v sin 0 cos ^ + v cos 0 cos <P) 

d v x 


* 

a 




in <P + cos 
°x 




(H-4) 


The partial derivatives of H- with respect to the R variables can be obtained 

in terms of the partial derivatives of H- with respect to the S variables by 

r 0 x y z 

using the Jacobian, J[(r,9, v , v ,v )/(x,y,z,v , \r ; v )]. The thirty-six 
terms of this Jacobian are given in Table III. They have been obtained by 
differentiating (II-2) with respect to the R variable and afterwards substituting 
the R variable by Expression (II-3). The final relations sought are: 


3p. . Q . m %*• cos 6 cos <P 0|x sin <P 1^,9 (n 

-rj- = -Tj- sin 0 cos Y + -os ; - tj f- n : n + (v cos 0 cos Y 

d x o r qd r °y r sin 0 r g r 


W 


<P //x . 1 , r a (0 *P cos 9 sin <p. 1 ^ , r . , n 

v T sin <P) - — ■-= — n (v cos 0 cos Y + v : zs ) + ^(v sin y + 

r V sin 0 r Q V 


0 CO s 0 

v sin <P — n) 

sin 0 


( II-5) 


fy. 9|x , 0 /n %x cos 0 sin <P cos <P 1 9^ , 9 tn 

^ sin 0 sin <P + w — + Wp r si~S + 7 77 (V sm *P cos 0 + 


<P 


1 % 


<P 


v cos <P) + n (- v cos 9 sm (P + v cos (P — 5 -) 

r a 0 sin 0 

d v 


9 X 1 ^ 


r7? 


9v^ 


m , r 0 cos 0 . 

cos Y (v 4* v — a) 

sin 0 
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3p. 3p. fl 3p. sin 0 1 ©P- © . 1 ©P- r 

J- - J- cos 9 - -to v sin 6 + n v sin 


^ = -g? 


■^0 


r r a r 


7"B 

dv 


sin 6 cos <P + -^75 cos © cos <jP - sin <P 

a v x a v r 0v e 


a. 


— = — - sin 0 sin 9 + — ^-n cos 0 sin 9 + — vn cos 9 

a v y a v r a v ° . ©v 


CO s 6 - * sin e 

a v z a v r ©v 0 


Substituting (II-5) into (II-4) and simplifying we obtain: 


(II-5) 


dr 

dt 


1 8[i 


2 (i a r 
p. d v 


d0 

dt 


e 

V 

rF* 


i apt 


2 * a v 


d^ 

dt 


1 , v^* _ 1 ©H- . 

in" 0 p 2 " * dj? 


r sm 


0 


dv r 1 ©p- . „0 f v v 1 ©P- 

^ 2 ~ rP- 

r p- a v 


<p . 


dt sln 6 [ - ? i/ 1 


1 ( v? 1 


0 r 

dv__ 1 j 1 
dt “ r [ P- 


Sp. 0 ,v r .1 9p. . 9 

*t’- v *2 ‘ * i7> + v r cos 


6 t-TTITB 


. 9 

dv 

dt " r sin 


> (■ 

in © [ P- 


■j© - sin 0 v<? (i - i -^) - r cos 0 v9 - ^7 

^ p. r a v 


v e 1 * l 
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Substituting the terms between parenthesis of the last three equations 
dr d0 d *P 

by — and — the set (II-6) will result in the ray tracing differential 
Equations (4.1). 


